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ABSTRACT 


The  statistical  covariance  technique  is  proposed  as  a  computer 
analysis  tool  for  missile  systems  error  propagation  studies.  It  is 
shown  that  on  representative  examples  the  statistical  covariance  tech¬ 
nique  performs  better,  based  on  speed  and  accuracy, than  the  traditional 
Monte  Carlo  approach  of  averaging  a  large  number  of  random  simulation 
runs.  Alt'  ough  exact  only  for  linear,  time-varying  systems,  the  covariance 
technique  bives  attractive  approximate  results  for  a  wide  range  of 
nonlinear  missile  systems.  The  extension  to  nonlinear  systems  utilizes 
incremental  equations  about  the  noise-free  solution.  Except  in  the 
extreme  cases  of  very  harsh  nonlinearities  and  high  input  noise  levels, 
the  statistical  covariance  technique  has  been  shown  to  yield,  for  the 
systems  considered,  results  comparable  to  1000  Monte  Carlo  runs.  More¬ 
over,  an  accuracy  prediction  procedure  has  been  developed  to  estimate 
the  error  to  be  expected  from  statistical  covariance  as  well  as  the 
number  of  Monte  Carlo  runs  required  for  comparable  accuracy.  Further¬ 
more,  an  automatic  sensitivity  program  is  suggested  for  utilizing  the 
statistical  covariance  technique  on  existing  simulations  without  repro¬ 
gramming.  It  is  shown  that  the  usual  number  of  Monte  Carlo  runs  yields 
very  inaccurate  results,  e.g.,  for  one  system  considered  25  runs  results 
in  an  error  of  30  percent  while  100  runs  gives  10  percent  error.  A 
total  computer  software  package  for  statistical  covariance,  including 
automatic  sensitivity,  accuracy  prediction,  and  miss  distance  programs, 
is  proposed  for  further  development. 
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1. 


Introduction 


a.  Noise  In  Missile  Systems 

The  consideration  of  noise  problems  is  an  inherent  part 
of  the  analysis  and  design  of  large-scale  missile  systems.  The  usual 
noise  problems  occur  when  either  system  parameters  or  disturbance  inputs 
vary  randomly.  Such  stochastic  variations  may  be  represented  as  random 
processes  or,  in  the  case  of  unknown  bias  effects,  as  random  variables. 
Random  disturbances  often  appear  as  noisy  signal  Inputs  or  as  measure¬ 
ment  noise.  Typical  noise  sources  in  missile  systems  are  radar  clutter, 
gyro  drift  and  bias,  boresight  and  alignment  errors,  and  wind  gust  dis¬ 
turbances.  In  all  cases,  random  variations  in  system  parameters  and 
external  inputs  result  in  errors  which  are  propagated  throughout  the 
missile  system. 

This  report  examines  the  problem  of  noise  propagation  in  large-scale  missile 
systems  from  two  basic  viewpoints .  The  first  of  these  is  the  Monte  Carlo  tech¬ 
nique,  which  exercises  a  computer  simulation  of  the  missile  system  with  a  large 
number  of  random  sample  functions .  The  second  approach  is  an  analytical 
method  referred  to  as  the  statistical  covariance  technique.  Numerical  results 
indicate  that  the  statistical  covariance  approach  is  an  effective  and  effi¬ 
cient  tool  for  most  missile  systems  analysis  situations  and,  inmost  practical 
instances,  a  significant  improvement  over  the  Monte  Carlo  method. 


b.  Background  on  the  Monte  Carlo  Method 

Previously,  the  primary  capability  for  Monte  Carlo  noise 
studies  was  based  on  the  use  of  analog  noise  generators.  Analog  sources 
of  noise  generally  depend  upon  a  gas  tube  in  a  magnetic  field.  Bekey 
and  Karplus  fl]  have  argued  that  the  specifications  of  analog  noise 
generators  are  usually  only  approximately  known  and  that  better  noise 
sources  are  available.  Other  noise  studies  rely  on  tables  of  random 
numbers  [2],  but  the  use  of  tables  requires  some  type  of  storage  device. 
More  recently,  digital  pseudo-random  number  generators  have  become  the 
standard  input  for  Monte  Carlo  simulations.  Chambers  [3],  Hull  and 
Dobell  [4],  MacLaren  and  Marsaglia  [5],  and  Gelder  [6]  are  among  those 
who  have  developed  mixed  congruential  and  multiplicative  recurrence 
formulas  for  generating  pseudo-random  numbers.*  These  formulas  yield 
number's  that  are  uniformly  distributed  on  the  unit  interval  (0,  1). 


*The  term  pseudo-random  is  appropriate  because  the  numbers  are  not 
actually  random  but  are  generated  by  deterministic  means.  This  deter¬ 
ministic  origination  has  the  inherent  advantage  of  providing  reproducible 
input  sequences  which  are  often  useful  for  computer  program  debugging. 
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The  problem  of  generating  digital  representations  of  correlated  con¬ 
tinuous  random  inputs  has  been  treated  in  a  report  by  Rowland  [7].  Tech¬ 
niques  for  designing  digital  shaping  filters  have  been  developed  in 
reports  by  Gujar  and  Kavanagh  [8],  Holliday  [9],  Baum  [10],  Broste  [11], 
and  Conner  [12].  Moreover,  digital  simulation  problems  encountered  when 
using  only  a  finite  sequence  of  random  data  are  discussed  in  reports 
by  Jansson  [13]  and  Brown  and  Rowland  [14].  A  comprehensive  study 
resulting  in  a  computer  software  package  for  Monte  Carlo  simulations  was 
performed  within  the  U.S.  Army  Missile  Command,  Redstone  Arsenal  by 
J.  L.  Harris  [15]. 

It  should  be  emphasized  that  any  statistical  information  related  to 
noise  propagation  in  missile  systems  simulations  must  agree  with  cor¬ 
rectly  implemented  Monte  Carlo  results. 


c .  Background  for  the  Statistical  Covariance  Technique 

The  statistical  covariance  approach  stems  from  the  use  of 
the  error  covariance  matrix  in  the  Kalman  filtering  equations  [16,  17]. 
Almost  every  new  book  in  the  last  five  years  on  stochastic  estimation 
and  control  includes  a  discussion  of  the  statistical  covariance  tech¬ 
nique  for  linear,  time-varying  systems  [18—23].  The  method  has  also 
been  applied  for  handling  variations  about  a  nominal  flight  path  obtained 
from  nonlinear  system  equations.  Kuhnel  and  Sage  [24]  applied  the 
covariance  technique  to  sensitivity  equations  about  the  nominal  flight 
path  caused  by  trajectory  initial  dispersions  and  parameter  variations. 
They  considered  a  thirty-third  order,  6  degree-of- freedom  homing  missile 
model  to  illustrate  the  usefulness  of  the  technique  to  a  realistic  situ¬ 
ation.  Irwin  and  Hung  [25]  applied  the  same  method  to  nonlinear  systems 
having  random  bias  inputs,  e.g.,  misalignment  or  drift,  errors.  They 
compared  the  direct  and  adjoint  approaches  for  evaluating  the  state 
covariance  matrix,  while  Kuhnel  and  Sage  used  only  the  adjoint  method. 
These  papers  demonstrated  that  the  covariance  technique  is  useful  for 
nonlinear  systems.  There  have  been  other  proposed  linearization  schemes 
which  yielded  varying  degrees  of  success  [26,  27].  These  approaches  are 
based  on  component  by  component  linearization  rather  than  the  more  power¬ 
ful  approach  of  linearizing  variational  equations  about  the  exact  nominal 
trajectory.  More  recently,  Brown  [28—30]  has  used  the  statistical 
covariance  approach  to  solve  trajectory  optimization  problems  for  non¬ 
linear  feedback  systems.  In  addition,  Clark  [31,  32]  has  developed 
related  results  using  improved  estimation  about  a  nominal  trajectory. 

The  statistical  covariance  technique  is  generally  recognized  as  the 
modern  approach  for  handling  error  propagation  in  guidance  and  control 
systems.  Its  application  to  missile  systems  analysis,  emphasizing 
advantages  and  limitations,  will  be  demonstrated  in  later  paragraphs  of 
this  report . 
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d. 


The  Error  Propagation  Problem 


The  analysis  problem  to  be  considered  here  may  be 
described  by  the  vector  differential  equation 


x  =  f (x,  r(t),  w(t),  t) 


(1) 


where 

x  =■  the  n-dimeusional  vector  of  system  variables 
r(t)  =  the  non-random  inputs 

w(t)  *  the  continuous  random  process  disturbance  inputs 

t  =  the  independent  variable  representing  time. 

For  simplicity  of  notation,  the  vector  r(t)  in  Equation  (1)  will  be 
suppressed  hereafter,  and  the  system  under  consideration  will  be 
expressed  as 


x  =  f(x,  w,  t)  .  (2) 

The  problem  is  to  determine  the  probability  density  function  of  the 
system  variables  x  as  a  function  of  time  t.  In  most  cases  of  interest, 
it  will  be  sufficient  to  find  only  the  mean  and  variance  of  x  (t)  . 


e .  Outline  of  the  Report 

Following  this  brief  introductory  background  material  and 
problem  definition,  Paragraph  2  describes  the  Monte  Carlo  approach  in 
detail  and  gives  numerical  results  for  a  simple  example.  It  is  shown 
that  a  very  large  number  of  simulation  runs  (up  to  1000)  are  needed  to 
yield  acceptable  accuracy.  The  statistical  covariance  technique  is 
introduced  in  Paragraph  3.  The  basic  approach  is  developed  for  linear 
systems  and  extended  for  an  approximate  analysis  of  nonlinear  systems. 
Paragraph  4  presents  program  refinements,  including  an  automatic  inter¬ 
rogation  program  for  determining  sensitivity  coefficients  in  existing 
large-scale  missile  simulations.  Paragraph  5  discusses  the  use  of  the 
statistical  covariance  technique  for  large-scale  systems  and  an  approxi¬ 
mate  sequential  algorithm  is  developed.  Finally,  Paragraph  6  summarizes 
the  conclusions  of  this  report  and  indicates  problem  areas  for  future 
consideration. 
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2.  Monts  Carlo  Studies 

The  traditional  approach  for  the  analysis  of  noise  propaga¬ 
tion  in  dynamical  systems  involves  Monte  Carlo-based  statistics.  Two 
elements  are  needed  in  the  digital  implementation  of  the  Monte  Carlo 
method.  The  first  of  these  is  the  digital  generation  of  a  sequence  of 
pseudo-random  numbers  for  input  to  the  given  system.  Several  multi¬ 
plicative  random  number  generators  are  investigated  in  this  paragraph. 

The  second  consideration  is  the  sampling  problem  -inherent  in  representing 
continuous  systems  and  signals  digitally.  The  modification  of  the  vari¬ 
ance  of  the  generated  random  sequence  as  a  function  of  the  integration 
step  size  is  necessary  to  handle  this  sampling  problem. 

This  paragraph  selects  a  particular  digital  random  number  generator 
for  use  in  subsequent  work  by  testing  several  popular  generators  on  a 
representative  system.  It  should  be  noted  that  even  this  "best"  genera¬ 
tor  yields  poor  results  when  only  100  to  200  Monte  Carlo  runs  are  used. 
Therefore,  the  case  for  implementing  the  statistical  covariance  tech¬ 
nique  to  be  described  in  Paragraph  3  is  strengthened. 


a.  Digital  Random  Number  Generators 

The  recurrence  formula  often  utilized  in  generating 
pseudo-random  numbers  has  the  form 

Zk+1  =  (Modulo  M)  ,  (3) 

which  is  referred  to  as  the  multiplicative  method.  The  scalar  constants 
A  and  M  may  be  selected  to  insure  good  statistical  properties.  Several 
requirements  for  selecting  A  and  M  are  given  in  a  report  by  Brown  and  Rowland  [14]  , 

20 

from  which  the  values  of  A  =  19971  and  M  =  2  are  highly  recommended 
for  use  with  a  starting  number  Zq  =  31571.  It  should  be  noted 

that  the  suggested  generator  yielded  pseudo-random  numbers  with  con¬ 
siderably  better  statistical  properties  than  several  other  widely  used 
generators,  including  the  Harris  generator  [15]. 

The  meaning  of  the  modulo  operation  in  Equation  (3)  is  that  the 
product  AZj^  is  first  divided  by  M  and  the  remainder  is  taken  as 

the  next  random  number.  Each  random  number  may  be  normalized  to  the 
unit  interval  by  dividing  by  M.  The  resulting  sequence  of  numbers  will 
be  approximately  uniformly  distributed  on  (0,  1) .  Moreover,  the  random 
number  sequence  satisfies  the  requirements  necessary  for  the  digital 
representation  of  Gaussian  white  noise. 


4 


Box  and  Muller  [33]  developed  a  direct  (or  exact)  approach  for 
transforming  two  independent  random  variables  from  the  same  uniform  dis¬ 
tribution  on  the  unit  interval  to  a  pair  of  independent  random  variables 
obeying  a  zero-mean,  unity-variance  normal  distribution.  The  uniform 
random  variables  Zt  and  Z.  are  related  to  the  normal  random  variables 
Yj  and  Y2  by 

1 

Y^  a  ^-2  10gfc  cos 

1 

Y2  =  ("2  loge  Zl)2  sin 

These  transformations  were  used  in  obtaining  some  of  the  numerical 
results  to  be  presented.  A  second  method  for  obtaining  normally  dis¬ 
tributed  random  numbers  is  to  sum  12  numbers  uniformly  distributed  on 
(0,  1)  and  subtract  the  sum  of  their  means.  This  approach,  referred  to 
as  the  central- limit  technique,  is  easier  to  apply  but  requires  a  larger 
number  of  random  numbers  for  averaging.  Both  methods  are  used  on  a 
second- order  system  later  in  the  following  section. 


2nZ2 

2*Z2  .  (4) 


b.  Numerical  Results 


Two  multiplicative  digital  pseudo-random  generators  v.ere 
tested  initially  on  the  XDS  Sigma  5.  Both  used  A  =  19971  with  a  starting 
number  Z  =  31571.  The  program  for  the  first  generator,  for  which 
20  ° 

M  =  2  =  1048576  as  suggested  by  Brown  [14],  is  shown  in  Figure  1.  Pro¬ 

gram  statements  4  through  7  provide  initialization  for  the  uniform  genera¬ 
tor  given  by  statements  13  through  21  and  28.  Statements  25  and  26 
transform  the  uniform  results  into  a  normal  distribution  according  to 
Equation  (4).  The  first  30  numbers  from  this  generator  are  given  in 
Figure  A-l  of  Appendix  A.  Unnormalized  integers  IX,  the  uniformly  dis¬ 
tributed  numbers  U,  and  the  normally  distributed  numbers  XNORM  are 
listed . 

*  31 

The  second  generator  differs  from  the  first  only  in  that  M  =  2  is 

used.  To  effect  this  change,  statements  14  through  21  of  the  program 

in  Figure  1  were  replaced  by  the  four  statements 

IX  =  IY 

IF(IX)  5,  5,  6 
5  IY  =  IY  +  2147483647+1 


6  U  »  IY*0 .4656612873E-9 


(5) 
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OPERATION  OF  A  MULTIPLICATIVE  DIGITAL  GENERATOR 


XMEAN=0 . 

SIG=1 . 

IX=31571 
DUM=0.1 
DO  99  K=l,30 

MULTIPLICATIVE  GENERATOR 

WITH  M=2  TO  THE  3 1ST  POWER 

IY=19971*IX 

IX=IY 

IF(IX)5,5,6 

5  IY=IY+2147483647+1 

6  U= IY *0. 46566 12873E-9 

TRANSFORMATION  TO  NORMAL 

Z=SQRT ( - 2 . 0*AL0G (DUM ) ) *S IG 
XNORM=Z*COS (6 . 28318*U)+XMEAN 

DUM=U 

PRINT  77,K,IX,U,XNORM 
77  FORMAT (6X,I4,8X, 112, 2F15. 6) 

99  CONTINUE 
STOP 
END 


Figure  1.  The  Multiplicative  Pseudo-Random  Number  Brown 
20 

(M  =  2  )  Generator 


Because  the  XDS  Sigma  5  is  a  32-bit  machine,  the  division  by  M  is  handled 

31 

automatically  when  M  is  selected  as  2  .  Therefore,  the  resulting  pro¬ 

gram,  as  evidenced  by  Equation  (5),  is  simpler.  The  latter  two  state¬ 
ments  in  Equation  (5)  represent  register  shifts  which  accomplish  the 
automatic  division.  Numerical  results  for  the  operation  of  the  second 
generator  are  given  in  Figures  A-2  and  A-3  of  Appendix  A. 

A  program  is  given  in  Figure  A-4  of  Appendix  A  for  examining  some 
of  the  statistical  properties  of  the  two  generators.  Sample  variance 


6 


and  means  are  plotted  in  Figures  2  and  3  as  functions  of  the  number  N  of 
pseudo-random  numbers  generated.  The  following  observations  may  be  made 
regarding  this  data: 

1)  The  prespecified  means  and  variances  are  obtained  as  N  becomes 
sufficiently  large. 

2)  For  finite  sequences  of  generated  numbers,  the  sample  means  and 
variances  vary  up  to  several  percent  about  the  prespecified 
values. 


N(THOUtANDS)  — » 

Figure  2.  Sample  Variances  for  Two  Multiplicative  Generators 


0123456789  10 

N  (THOUSANDS)  . . .  — . 

Figure  3.  Sample  Means  for  Two  Multiplicative  Genrators 

Another  important  statistical  property  to  be  considered  in  random 
number  generation  is  the  autocorrelation  function.  It  is  noted  in  a 
report  by  Brown  and  Rowland  [14]  that  the  serial  autocorrelation  coeffi¬ 
cients  for  independent  Gaussian  random  sequences  should  be  asymptotically 
normally  distributed  with  mean  zero  and  variance  1/N.  The  Brown  raulti- 
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plicative  generator  with  M  =  2  satisfies  the  hypothesis  of  independence 
with  a  high  degree  of  probability  [14], 
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c .  Discrete  Representations  for  Continuous  Noise  Processes 

The  use  of  random  numbers  as  disturbance  input  sample 
functions  for  dynamical  systems  requires  special  considerations.  It  is 
necessary  to  determine  discrete  representations  for  continuous  noise 
processes  [7] .  If  pseudo-random  numbers  are  held  constant  over  one 
sampling  period  H,  then  the  corresponding  autocorrelation  function  is  the 
triangular  form  shown  in  Figure  4.  It  is  required  that  this  triangular 
function  represent  as  closely  as  possible  the  impulse  function  for  the 
continuous  case.  Equating  the  power  spectral  densities,  i.e.,  Fourier 
transforms  of  the  autocorrelation  functions,  of  discrete  and  continuous 
cases,  yields 


Qw 

Qwd  -  T  ’  (6> 

where  and  are  the  variances  of  the  discrete  and  continuous  cases, 

respectively.  Together  with  the  proper  utilization  of  a  multiplicative 
pseudo-random  number  generator,  the  relationship  in  Equation  (6)  is  the 
most  important  element  in  Monte  Carlo  simulations. 


AUTOCORRELATION 

FUNCTIONS 


POWER  SPECTRAL 
DENSITIES 


CONTINUOUS 
WHITE  NOISE 


DISCRETE 
WHITE  NOISE 


Figure  4.  Discretization  Effects  for  White  Noise 


d .  Simulation  Results 

Consider  the  t ime- invariant ,  second-order,  linear  system 

described  by 


x 


1 


=  x 


2 


*2  ■  -2x^  -  3x2  +  w 


(7) 
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with  x, (0)  =  x2(0)  =  0.  Equation  (7)  gives,  in  phase  variable  form, 
the  state  equations  of  the  system  with  transfer  function 


H(s) 


_ 1 _ 

(s  +  l)(s  +  2) 


(8) 


where  x^  represents  the  output  x  and  x2  represents  x.  The  input  w  is  a 

zero-mean,  unity -variance,  gaussian  white  noise  signal  applied  for  all 
t  ^  0.  The  problem  is  to  investigate  the  performance  in  Monte  Carlo 
simulations  for  Equation  (7)  by  using  several  multiplicative  pseudo¬ 
random  number  generators . 
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The  computer  program  used  for  testing  the  Brown  (M  =  2  )  generator 

with  the  exact  nonlinear  transformation  in  Equation  (4)  is  included  as 
Figure  A-5  of  Appendix  A.  The  program  shown  was  changed  (statement  10) 
to  allow  25,  50,  100,  200,  500,  and  1000  simulation  runs  to  be  ensemble- 
averaged  at  half-second  intervals  from  t  =  0  to  t  =  5  seconds.  The 
variance  of  the  output  x  is  plotted  as  a  function  of  time  t  in  Figure  5 

20 

for  the  Brown  (M  =  2  )  generator  using  the  exact  transformation  in 

Equation  (4) .  The  exact  variance  solution  is  determined  by  convolution 
in  Appendix  8*.  Large  errors  were  obtained  for  all  curves  except  the 
case  using  1000  runs.  ' 

Figures  6  and  7  show  the  average  percent  error  obtained  for  the 
output  variance  of  Equation  (7)  as  a  function  of  the  number  of  runs  for 
several  generators.  Using  the  basic  program  in  Figure  A-5  in  Appendix 
A,  statements  27  through  34  were  modified  according  to  Equation  (5)  to 
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give  corresponding  data  for  the  Brown  (M  =  2  )  generator.  Two  new 

multiplicative  generators  were  also  testea.  One  of  these,  the  SAM-D 
generator,  has  been  used  in  guidance  simulation  studies  for  a  large-scale 
air  defense  missile  system.  The  SAM-D  generator  has  a  multiplier  A  of 
899  and  a  starting  number  Zq  =  1073741823.  The  other  generator  [34], 

which  has  been  widely  used  on  the  IBM  360,  has  a  multiplier  A  and  starting 
number  Zq  =  1220703125.  Figure  6  shows  results  for  these  four  generators 

obtained  by  using  the  exact  nonlinear  transformation  in  Equation  (4)  from 
uniform  to  normal.  Figure  7  shows  results  using  the  central-limit  approach 
approach  of  summing  12  numbers  from  the  uniform  generators.  It  should  be 

20 

noted  that,  while  the  Brown  (te=2  )  generator  had  an  average  error  of  only 

1.5  percent  for  1000  runs,  the  largest  error  obtained  at  the  10  equally- 
spaced  test  points  was  2.7  percent.  Therefore,  results  for  even  1000  runs 
may  be  considered  to  be  only  marginally  acceptable  in  some  applications. 


The  statistical  covariance  technique  to  be  described  in  Paragraph  3  also 
yields  the  exact  solution.  Moreover,  the  statistical  covariance  technique 
is  much  easier  to  apply  as  well  as  being  more  amenable  for  machine  com¬ 
putation  than  the  convolution  approach. 
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OUTPUT  VARIANCE 


Figure  5.  Monte  Carlo  Results  for  the  System  in  Equation  (7)  Using 
20 

the  Brown  (M  =  2  )  Generator  with  the  Exact  Nonlinear 

Transformation  to  Normal 
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AVERAGE  PERCENT  ERROR  ON  OUTPUT  VARIANCE 


NUMBER  OF  RUNS  <LOG10  SCALE)  - - 

Figure  6.  Average  Percent  Error  on  Output  Variance  Versus  the  Number 
of  Monte  Carlo  Runs  for  Generators  Using  the  Exact  Nonlinear 
Transformation  in  Equation  (4) 
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Step  Sizes  Versus  Number  of  Runs 


Consideration  should  also  be  given  to  the  problem  of  how 
to  use  a  fixed  number  of  random  numbers  in  a  Monte  Carlo  simulation. 

For  example,  the  question  arises  as  to  whether  the  error  in  the  output 
variance  would  be  smaller  by  using  a  large  step  size  with  many  runs  or 

a  small  step  size  with  fewer  runs.  If  N  is  the  number  of  sample  points 

P 

per  run  and  ^  is  the  number  of  runs  to  be  ensemble-averaged,  then  the 
error  in  the  output  variance  obviously  decreases  as  either  or 
increases .  Note  that  the  produce  N^N^  is  the  total  number  of  random 
numbers  to  be  used.  The  problem  is  to  determine  the  values  of  and  N^, 
when  the  product  N^N^  is  fixed,  to  yield  the  smallest  possible  error. 


Figure  8  shows  plots  of  average  error  for  the  output  variance  as  a 
function  of  the  step  size  H  for  the  system  in  Equation  (7)  when  the  first 
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100,000  random  numbers  are  used  from  the  Brown  (M  =  2  )  generator.  For 

example,  the  value  at  H  =  0.01  for  the  5-second  problem  was  obtained  by 

using  N  =  500  and  N  =  200.  Two  methods  of  averaging  the  error  vari- 
P  ^ 

ance  data  at  half-second  intervals  were  used.  The  first  treated  all  data 
points  equally,  but  the  second  approach  weighted  the  data  according 
to  the  number  of  random  numbers  which  had  been  put  into  the  system 
prior  to  that  particular  problem  time.  For  example,  the  data  at 
t  =  3.0  seconds  was  weighted  with  a  factor  of  only  60  percent  as  large 
as  the  data  at  t  =  5.0  seconds,  because  only  60,000  random  numbers 
were  applied  to  the  system  up  to  t  =  3.0  seconds  while  100,000  were 
used  on  the  range  t  =  0  to  5.0  seconds.  The  curve  in  Figure  7  shows 
that  as  large  a  step  size  as  possible  should  be  used  without  violating 
conditions  for  which  truncation  errors  become  significant  in  integrating 
the  system  equations.  In  this  example,  the  maximum  step  size  under 
those  conditions  was  H  =  0.5  second.  Note  that  the  use  of  a  somewhat 
smaller  step  size  gave  much  larger  average  errors  in  the  output 
variance . 


f .  Sumnary 


A  detailed  description  has  been  given  on  the  digital 
implementation  of  the  Monte  Carlo  method.  Discretization  problems  and 
pseudo-random  number  generation  have  been  discussed,  and  several  multi¬ 
plicative  generators  have  been  tested  on  a  second-order  linear  system. 

20 

Simulations  showed  that  the  Brown  (M  =  2  )  generator  performed  consis¬ 

tently  well  in  a  variety  of  tests.  An  interesting  result  is  that  for  max¬ 
imum  efficiency  in  Monte  Carlo  simulations  the  largest  step  size  compati¬ 
ble  with  integration  accuracy  should  be  used  with  a  large  number  of  com¬ 
puter  runs  for  ensemble-averaging.  Inaccuracies  of  the  Monte  Carlo 
method,  as  well  as  the  large  amount  of  computer  time  required,  point  to 
the  need  for  the  statistical  covariance  technique  of  Paragraph  3. 
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Figura  8.  Determining  the  Optimum  Utilization  of  100,000  Random  Numbers 
i:  Monte  Carlo  Simulations  for  System  in  Equation  (7) 

3.  The  Statistical  Covariance  Technique 

An  analytical  method  is  developed  in  this  paragraph  and  com¬ 
pared  on  a  nonlinear  system  with  the  Monte  Carlo  approach  described  in 
Paragraph  2.  The  statistical  covariance  technique  is  derived  as  an  exact 
solution  for  linear,  time-varying  systems  and  then  extended  for  an  approx¬ 
imate  analysis  of  the  error  propagation  problem  in  nonlinear  systems. 

The  method  is  also  applied  to  the  problems  of  miss  distance  studies  and 
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error  budgeting.  The  primary  conclusion  of  this  paragraph  is  that  the 
statistical  covariance  technique  is  valid  for  a  wide  range  of  input 
noise  signals  and  parameters  in  nonlinear  system  analysis. 


a.  Derivation  of  Statistical  Covariance  Results 


Let  the  system  under  consideration  be  described  by  a 
linear,  time-varying  vector  differential  equation  of  the  form 

x  =  Ax  +  Bw  ,  (9) 


where 


x  =  the  n-dinensional  state  vector 
w  =  an  m-dimensional  input  vector 
A  =  an  n  by  n  matrix 
B  =  an  n  by  m  matrix 

The  input  vector  w  has  elements  that  are  zero-mean  white  noise  with 
a  covariance  matrix  given  by 

Ejw(t)  wT(t)J  A  QwS(t  -  T)  •  (10) 

The  matrix  Q^,  as  well  as  the  matrices  A  and  B,  may  be  time-varying. 

Let  the  covariance  matrix  of  the  state  x  be  defined  by 

p(t)  A  e|x  (t)  xT(t)|  ,  (11) 

where  T  denotes  the  transpose. 

The  problem  is  to  determine  P(t)  in  terms  of  A,  B,  and  Q  .  Three 

approaches  to  the  problem  will  be  used.  The  first  two  apply  directly 

to  the  continuous  system  in  Equation  (9)  and  the  third  pertains  to  a 

discretized  version  of  the  system. 

(1)  The  Integral  Solution.  The  first  derivation  of  an 
expression  for  P(t)  utilizes  the  integral  solution  for  x(t).  From 
Equation  (9), 
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Equation  (15)  is  the  desired  integral  solution  for  P(t).  However,  by 
forming  P(t)  from  Equation  (15)  and  using  the  relationship 

=  AJ(t,  r)  ,  (16) 

the  result  may  be  expressed  as  the  following  matrix  differential 
equation: 

P  =  AP  +  PAT  +  BC^B1  .  (17) 


This  differential  form  of  the  statistical  covariance  equation  is  the 
primary  result  of  Paragraph  2. 

(2)  The  Differential  Solution  Directly.  Using  Equation 

(9)  directly, 

P  =  ^  j^|x(t)  xT(t)(]  =  e|xxT  |  +  e|xxT  | 

=  E  |(ax  +  Bw)  xT)  +  e|x(Ax  +  Bw)T| 

=  A  E jxxT  |  +  B  e|wxT  |  +  e|xxT  J  Ax  +  e|xwTJbT 

=  AP  +  PAT  +  B  e{wxT J  +  e|xwT|  BT  .  (18) 

Substituting  Equation  (12)  into  Equation  (18)  for  x(t)  and  again  applying 
the  sifting  property  of  the  delta  function  to  perform  the  integrations, 
Equation  (17)  is  obtained  [18]. 

(3)  The  Discrete  Version.  Equation  (12)  may  be  written 
in  the  discrete  version  as 


~k-f  1  *  \  *lc  +  V^k 


(19) 


Therefore, 
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which  is  the  discrete  equivalent*  of  Equation  (17). 


(20) 


b.  Approximate  Analysis  of  Nonlinear  Systems 

The  application  of  the  statistical  covariance  technique 
to  nonlinear  systems  can  be  achieved  as  an  approximate  analysis.  Con¬ 
sider  small  variations  Sx(t),  caused  by  the  input  disturbance  noise 
w(t),  about  the  (noise-free)  nominal  trajectory  xN(t),  as  shown  in 

Figure  9.  These  5x  variations  are  assumed  to  be  sufficiently  small  to 
satisfy  linear  perturbation  equations,  i.e., 

5x(t)  =  x(t)  -  xN(t)  (21) 


and 


6x  =  A(t)5x  +  B(t)  w(t) 


(22) 


where 


dl 

-  k 


B(t)  = 


df 

dw 


x  =  %(*> 
w  =  \ 


=  xN<t) 

w  =r 


(23) 


S  =  Tw 


♦Recall  that  and  ,  the  covariance  matrices  for  the  continu¬ 
ous  and  discrete  cases,  respectively,  are  related  by  =  0^/H,  where 
H  is  the  discretization  interval  as  given  in  Equation  (6)  of  Paragraph  2. 
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Figure  9.  A  Sketch  Showing  Incremental  Variations  &x(t) 
About  a  Nominal  Trajectory  xN(t) 


denotes  the  mean  of  w, (which  has  been  assumed  zero  in  the  problem 

formulation) .  The  nonlinear  vector  function  f  is  given  in  Equation  (2) . 
The  only  approximation  made  is  that  the  second  and  higher  order  terms 
in  6x  are  negligible  when  compared  with  the  linear  terms  which  appear 
in  Equation  (22).  This  approximation  is  valid  when  the  5x  variations  are 
sufficiently  small. 

(1)  Numerical  Results  for  Nonlinear  Systems.  The  fol¬ 
lowing  numerical  results  are  given  to  provide  insight  into  the  practi¬ 
cality  of  using  the  statistical  covariance  technique  for  nonlinear 
systems.  In  particular,  the  ranges  of  noise  signals  and  parameters  of 
the  system  nonlinearity  are  determined  for  which  the  statistical  covari¬ 
ance  technique  is  recommended.  Comparisons  are  made  with  25,  50,  100, 
200,  and  1000  Monte  Carlo  runs. 


(a)  Example  1  —Consider  the  second-order  nonlinear 
system  shown  In  Figure  10.  The  system  equations  are 

2 

Xj^  =  -  2x^  +  x2  +  a  x2  sign  ^x^ 

x2  =  -x2  +  w(t)  ,  (24) 

where  w(t)  Is  a  Gaussian  white  noise  process  applied  for  t  ^  0.  The 
process  w(t)  has  zero  mean  and  variance  Q^.  The  initial  condition  c*n 

x^  is  zero,  but  x^(0)  is  allowed  to  vary  in  parts  of  this  first  example. 

The  statistical  covariance  technique  and  the  Monte  Carlo  method 
were  progratmted  for  system  Equation  (24)  in  the  single  computer  program 
given  in  Figure  C-l  of  Appendix  C.  The  step  size  of  H  =  0.05  was 
selected  according  to  the  pseudo-random  number  utilization  criterion 
developed  in  Paragraph  2 .  Note  that  for  a  =  0  and  x2(0)  =  0  the  system 

in  Figure  10  is  identical  to  system  Equations  (7)  and  (8)  of  Paragraph  2, 
except  that  a  different  set  of  state  variables  has  been  chosen  in 
Equation  (24) .  Figure  11  shows  a  sketch  of  typical  Monte  Carlo  runs 
along  with  the  one-sigma  contours  from  the  statistical  covariance  results 
for  a  =  0.05,  x2(0)  =  0.1,  and  =  0.05.  As  the  number  of  Monte  Carlo 

runs  becomes  large,  i.e.,  on  the  order  of  1000  runs,  the  one-sigma  con¬ 
tours  from  the  statistical  covariance  technique  are  approximately  the 
same  (within  2  percent)  as  the  one-sigma  contours  computed  from  the 
sample  functions. 

The  validity  of  the  statistical  covariance  technique  for  several 
values  of  a  and  is  examined  in  Figures  12  and  13.  In  each  case,  the 

statistical  covariance  result  agreed  within  2  percent  with  1000  Monte 
Carlo  runs  for  small  values  of  the  parameters.  As  a  and  were  increased, 

the  error  in  the  statistical  covariance  solution  also  increased.  These 
curves  are  used  later  to  estimate  the  accuracy  expected  from  the  statis¬ 
tical  covariance  technique  by  examining  the  nonlinear  system  equations  directly  . 

Comparisons  with  a  limited  number  of  Monte  Carlo  runs  are  given  in 
Figures  14  and  15.  The  first  of  these  figures  compares  time  plots  of 
the  output  variance  for  25,  50,  100,  200,  and  1000  runs  with  the  sta¬ 
tistical  covariance  results.  Figure  15  shows  averaged  results  as  percent 
error  on  the  output  variance  versus  the  input  noise  covariance  C^.  All 

error  curves  are  based  on  1000  runs  as  the  standard  for  comparison.  In  view 
of  the  results  of  Figures  6  and  7  in  Paragraph  2,  it  should  not  be  too 
surprising  to  find  that  even  200  runs  yield  errors  of  approximately  8 
percent,  while  only  25  runs  yield  approximately  25  percent  error. 
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SYSTEM 

NONLINEARITY 


Figure  10.  The  Open- Loop  System  and  System  Nonlinearity 
Described  by  Equation  (24) 


An  interesting  problem  was  encountered  when  varying  X£(0)  for  this 

example.  It  was  found  that,  for  a  =  0.05  and  =  0.1,  the  percent 

difference  between  the  statistical  covariance  result  and  1000  Monte  Carlo 
runs  increased  rapidly  as  X£(0)  was  increased.  However,  for  reasons 

explained  later,  it  was  suspected  that  the  1000  Monte  Carlo  simulation 
runs  for  large  x,(0)  were  not  giving  correct  results.  To  prove  this 


Figure  12.  Variations  in  Statistical  Covariance  Results  as  a  Function 
of  the  Amount  of  System  Nonlinearity 

conjecture,  the  linear  system  (a  =  0)  was  investigated  by  varying  X2(0) 

as  shown  in  Figure  16.  After  rechecking  the  computations  to  verify  that 
step  size  and  integration  accuracy  were  acceptable,  it  was  concluded  that 
the  inaccuracy  of  the  digital  random  number  generator  must  be  the  cause 
of  the  errors.  In  particular,  the  non-zero  means  and  nonunity  variance 
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TIME  (mc) 


Figure  16.  Time  Plots  of  the  Output  Variance  for  Various  x2<0) 
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Xj  =  -  2Xj^  +  x2  +  ax 2  Sign  ^x^ 

x2  =  -Xj  -  x2  +  w(t)  .  (25) 

Note  that  Equation  (25)  is  the  system  Equation  (24)  of  Figure  10  with  the 
addition  of  a  negative  unity  feedback  loop.  Results  similar  to  those  in 
Figure  15  for  the  open-loop  case  were  obtained  for  Equation  (25).  These 
curves  are  shown  in  Figure  18. 


Figure  17.  Closed-Loop  System  Description  for  Equation  (25) 


(2)  Predicting  Statistical  Covariance  Accuracy.  It  would 
be  desirable  to  be  aole  to  predict  in  advance  the  accuracy  of  the  sta¬ 
tistical  covariance  technique  for  nonlinear  systems.  An  exact  prediction 
of  the  expected  accuracy  is  not  possible  because  no  exact  analytical 
solution  can  be  found,  in  general,  for  the  output  variance  of  nonlinear 
systems.  However,  the  result  from  a  large  number  of  Monte  Carlo  runs 
may  be  regarded  as  "exact"  for  the  purpose  of  accuracy  prediction,  but 
even  then  (as  shown  in  Paragraph  2)  some  inaccuracy  is  present.  The 
reason  for  using  the  statistical  covariance  technique  is  to  avoid  the 
time-consuming  Monte  Carlo  approach. 

Suppose  the  Monte  Carlo  runs  had  been  made  for  one  particular 
design  condition  (parameter  setting)  of  a  given  missile  system.  Using 
this  information,  the  following  procedure  could  be  used  to  estimate  the 
accuracy  of  the  statistical  covariance  technique  for  sufficiently  small 
changes  in  the  parameter  settings.  As  a  particular  example  to  illustrate 
the  procedure,  consider  the  exact  incremental  equation  associated  with 
Equation  (24),  i.e., 
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Figure  18.  Monte  Carlo  and  Statistical  Covariance  Results  Versus 
for  the  Closed-Loop  System  for  Equation  (25) 
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&x  +  2a8  x 
N  1 

&x2  =  - &x2  +  w(t)  .  (26) 

Suppose  that  the  nonlinear  term  in  Equation  (26)  is  required  to  be  not 
greater  than  k  percent  of  the  corresponding  linear  terms,  i.e., 


6x, 


=  -  28xt  + 


[i  +  2a|x2j 
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Squaring  Equation  (27)  and  taking  expected  values  yields 

*4.,  *['  +  2oN]  V, 

N  Z 


+  4 

1  +  2  a|x2j 

]  |e/sx  &x  l 

- 

J  N 1  l  I  | 

(27) 


(28) 
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Note  that  E<8x0  >  has  been  approximated  by  3  a  ,  which  is  exact  in  this 
V  *■ )  ®X2 

case  because  8x2  is  Gaussian.  Using  the  steady-state  values  of  the  vari¬ 
ance  terms  obtained  from  the  linear  case  (a  =  0)  yields 


8Xj  12 


12 


j8x2  -  2 


E  |&x^  6x2 


}  =  T 


(29) 


Substituting  Equation  (29)  into  Equation  (28)  gives,  after  simplifications, 


/  \2  2 
I  +  !  |l  *  2ox2(t)|  +|  |l  +  2o»2(t)| 

\  /  N  N 


(30) 


The  equality  in  Equation  (30)  is  plotted  in  Figure  19,  which  shows  that 
as  either  a  or  increases,  the  percent  k  of  the  first  incremental 

equation  in  Equation  (26)  caused  by  the  nonlinear  term  increases  rapidly. 
Using  the  information  in  Figure  19  together  with  Figures  12  and  13,  the 
percent  error  in  the  output  variance  as  a  function  of  the  parameter  k 
may  be  plotted.  The  sketch  for  varying  a  and  is  shown  in  Figure  20. 

If  k  is  less  than  5  percent,  then  the  error  in  the  output  variance  is  also 
less  than  5  percent.  However,  for  k  =  20  percent,  the  error  in  the  out¬ 
put  variance  is  approximately  30  percent.  If  a  and  are  such  that  k  is 
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Figure  19.  Plots  of  k  Versus  a  and  Q 

w 


approximately  20  percent,  then  the  statistical  covariance  technique  com¬ 
pares  in  accuracy  to  approximately  25  Monte  Carlo  runs  (30  percent  error) . 
However,  if  k  =  5  percent,  then  the  accuracy  of  the  statistical  covariance 
technique  is  better  than  200  Monte  Carlo  runs.  Therefore,  k  may  be  com¬ 
puted  in  advance  from  the  incremental  equations  to  determine  the  expected 
accuracy  and  the  number  of  Monte  Carlo  runs  which  would  yield  approximately 
the  same  accuracy  as  the  statistical  covariance  technique  for  the  given 
nonlinear  system. 


k - ► 

Figure  20.  Plots  of  Percent  Error  for  Statistical  Covariance  Versus  k 


These  observations  on  the  accuracy  of  the  statistical  covariance 
technique  as  a  function  of  the  quantity  k  are  correct  only  for  the 
single  example  previously  considered.  However,  it  can  be  expected  that 
other  second-order  systems  with  parameters  sufficiently  near  those  of  the 
previous  example  would  yield  results  with  corresponding  accuracy.  In 
particular,  it  should  be  expected,  as  shown  in  Figure  20,  the  average 
percent  error  in  output  variance  would  be  on  the  order  of  1.5  times  the 
value  of  k.  Moreover,  some  useful  information  would  be  obtained  even  if 
this  error  varied  by  as  much  as  1  to  2  times  k.  However,  variations  of 
from  50  to  100  times  k  would  be  unexpected. 
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Monte  Carlo  simulation  experience  is  usually  available  on  those 
missile  systems  where  noise  disturbances  have  been  a  problem.  Curves 
similar  to  those  in  Figure  20  can  be  plotted  for  the  particular  missile 
system  being  considered.  As  stated  previously,  these  curves  can  be  used 
to  yield  approximate  estimates  of  the  accuracy  of  the  statistical 
covariance  technique  in  given  situations. 


c.  Miss  Distance  Studies 


It  is  important  to  be  able  to  apply  the  results  of  the 
statistical  covariance  technique  to  miss  distance  studies  in  missile 
systems  analysis.  The  only  approximation  made  in  extending  the  basic 
covariance  method  to  nonlinear  systems  is  that  the  incremental  variations 
about  the  nominal  flight  conditions  obey  a  linearized  differential  equa¬ 
tion.  As  previously  stated,  this  approximation  is  valid  for  sufficiently 
small  noise  disturbances  or  small  nonlinear  variations.  However,  to 
apply  the  statistical  covariance  technique  to  the  missile  miss  distance 
problem,  one  further  assumption  is  needed.  This  assumption  is  that  the 
incremental  system  variables  &x  are  Gauss ianly  distributed  with  a  mean 
of  zero  and  a  covariance  given  by  P(t)  from  the  statistical  covariance 
solution. 

There  are  several  reasons  to  believe  that  the  Gaussian  assumption 
is  a  good  one.  First,  recall  that  many  of  the  noise  input  disturbances 
encountered  in  missile  systems  are  approximately  Gaussian.  While  Gaus¬ 
sian  inputs  to  nonlinear  systems  do  not  result  precisely  in  Gaussianly 
distributed  system  states,  the  deviation  from  the  Gaussain  condition  is 
caused  only  by  the  inexactness  of  the  approximation  in  linearizing  the 
perturbation  equations  in  5x.  If  these  linearized  equations  were  exact, 
then  6x  would  be  Gaussian  for  Gaussian  inputs.  A  second  reason  for 
giving  credence  to  the  Gaussian  assumption  is  the  Central-Limit  Theorem 
effect,  which  states  that  as  the  number  of  independent  inputs  becomes 
larger,  the  total  effect  on  the  system  tends  toward  a  Gaussian  condition. 
Finally,  experience  has  shown  that  the  convolution  property  of  linear 
subsystem  equations  inevitably  makes  the  subsystem  outputs  much  closer 
to  being  Gaussianly  distributed  than  whatever  input  distribution  was 
applied. 

Applying  the  Gaussian  assumption,  the  joint  probability  density 
function  for  the  incremental  state  variation  &x  may  be  written  as 


p(Sx)  =  - ^ - - -  exp  j^-  -r  SxT  P  1  Sx  J.  (31) 

(2n)n  ydet  P(t) 
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If  the  statistical  covariance  result  P(t)  is  used  only  to  compute 
the  miss  distance  m^,  not  all  components  of  P(t)  are  needed.  In  that 

case,  a  submatrix  is  formed  by  selecting  from  P(t)  those  components 
included  in 


(32) 


where 


$ E 


In  Equation  (32),  6y  and  6z  represent  the  incremental  variations  in 
orthogonal  components  of  the  miss  distance  (Figure  21).  Therefore,  if 
Sy  and  &z  are  assumed  to  be  jointly  normal  (Gaussian),  the  resulting 
joint  density  function  is 


If  Equation  (31)  holds,  then  so  must  Equation  (33).  However,  Equation 
(33)  may  be  true  without  Equation  (31)  being  true,  i.c.,  in  miss  distance 
studies  using  the  statistical  covariance  method,  it  is  necessary  to 
assume  only  that  5y  and  5z  are  jointly  normal;  the  remaining  system 
states  may  not  be  Gaussianly  distributed . 

In  calculating  the  circular  error  probable  (CEP)  in  miss  distance 
studies,  the  probability  density  function  of  the  miss  distance,  md,  must 

be  determined.  The  random  variable  m^  may  be  expressed  either  as 
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2  2  1 
where  o^  1  and  o^l  are  the  variances  of  the  transformed  variables  &y  and 

5z\  respectively.  Moreover,  the  probability  distribution  function 
P(inj)  is  given  by 


a 

P(md)  =  /  P(”d)  d”d 


and  the  CEP  is  defined  as  that  value  of  m^  for  which 


’(■“)  k  ■  * 


It  should  be  pointed  out  that  the  definite  integral  in  Equation  (37) 
can  be  solved  analytically  only  for  the  case  in  which 


W  -  V1  • 

Otherwise,  the  integration  may  be  performed  numerically  on  the  digital 
computer.  If 

"sV  -  V'  - 


then  p(m^)  In  Equation  (36)  is  the  Rayleigh  probability  density  function 
given  by 


md 

"  md2 

2  exp 

„  2 

.  &yl 

0 

for  m,  5  0 
d 


otherwise 


Moreover,  Equation  (38)  yields 


'("«>)  -  1 


from  which 


p("a)  ■  I 

when 

”d  -  °6yl  72  ln  2  -  %yl 

Therefore,  the  CEP  has  been  obtained  analytically  for  this  special  case 
as  1.17  cr  i.  The  probability  density  function  p(m,)  is  shown  in  Figure 

oy  2  2  a 

22.  Whenever  *  cr^i,  resu^t^,n8  probability  density  curve 

obtained  via  numerical  integration  has  a  similar  shape,  though  not  pre¬ 
cisely  a  Rayleigh  curve. 


Figure  22.  The  Rayleigh  Probability  Density  Function 


In  summary,  results  from  the  statistical  covariance  technique  can 
be  used  to  determine  statistical  information  about  missile  systems 
miss  distance.  The  procedure  is  to  transform  the  correlated  orthogonal 
components  of  the  miss  distance  into  an  uncorrelated  set  of  orthogonal 
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components  &y*  and  6z*.  The  density  and  distribution  functions  of 

are  found  from  Equations  (36)  and  (37)  in  terms  of  <r  l  and  a  i.  Sub- 

5y  gz 

sequent ly,  the  CEP  can  be  determined  by  using  the  distribution  function. 


d.  Error  Budgeting 

The  statistical  covariance  technique  can  be  applied  to  the 
problem  of  error  budgeting  for  determining  what  part  of  the  output 
variance  is  caused  by  any  specific  noise  source.  For  example,  if  the 
noise  source  is  the  result  of  errors  in  a  certain  measuring  device,  such 
as  a  particular  sensor,  then  it  would  be  advantageous  to  know  its  effect 
on  the  missile  miss  distance.  The  term  "error  budgeting"  refers  to  the 
synthesis  problem;  the  allowable  miss  distance  or  "error"  is  given  and 
the  objective  is  to  budget  this  allowable  error  among  several  contribu¬ 
ting  sources. 

The  statistical  covariance  technique  uses  linearized  variational 
(or  incremental)  equations  about  a  nominal  flight  path.  As  a  result, 
the  superposition  principle  can  be  applied  using  each  independent  input 
noise  source  individually  with  all  others  set  to  zero.  Appropriate 
allocations  can  be  made  on  this  basis  by  adjusting  weighting  coefficients 
on  the  input  noise  sources  until  an  acceptable  combination  is  determined. 


e.  Summary 


The  statistical  covariance  technique  has  been  developed 
in  the  paragraph  for  linear,  time-varying  systems.  The  practical  use 
of  the  technique  in  an  approximate  analysis  of  nonlinear  systems  has 
been  demonstrated  on  a  second-order  example.  From  the  exact  nonlinear 
incremental  equations ,  the  accuracy  expected  to  be  obtained  by  using  the 
statistical  covariance  technique  may  be  predicted  for  similar  systems. 
Furthermore,  in  such  cases,  an  estimate  of  the  number  of  Monte  Carlo 
runs  needed  to  achieve  this  same  accuracy  may  be  determined. 

In  stannary,  the  most  dramatic  result  of  Paragraph  3  is  that  not 
only  does  the  statistical  covariance  technique  yield  useful  results  for 
nonlinear  systems  but  an  estimate  of  its  accuracy  and  the  number  of 
Monte  Carlo  runs  needed  for  comparable  accuracy.  In  Paragraph  4,  pro¬ 
gram  refinements,  such  as  an  automatic  program  for  calculating  incre¬ 
mental  equations,  are  described  in  detail.  Subsequently,  in  Paragraph  5, 
a  sequential  algorithm  is  developed  to  apply  statistical  covariance  more 
efficiently  to  large-scale  missile  systems. 
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4.  Program  Refinements  for  the  Statistical  Covariance  Technique 

a.  Introduction 

The  use  of  the  basic  statistical  covariance  technique  on 
simple,  low-order,  linear  systems  is  straightforward.  The  method  dis¬ 
cussed  in  Paragraph  3  is  much  more  satisfactory  than  the  alternate 
approach  (Monte  Carlo)  of  ensemble-averaging  a  large  number  of  simula¬ 
tion  runs.  Moreover,  statistical  covariance  applies  to  a  wide  range  of 
nonlinear  problems,  as  described  in  Paragraph  3.  If  the  combination  of 
system  nonlinearities  and  noise  levels  remain  within  acceptable  limits, 
then  the  results  of  the  statistical  covariance  technique  are  important 
in  missile  system  analysis  and  design  studies. 

This  paragraph  examines  the  overall  usefulness  of  the  statistical 
covariance  technqiue  as  a  computer  analysis  tool  for  large-scale  missile 
problems.  Initially,  the  basic  objectives  are  considered  with  special 
emphasis  on  their  interdependence.  Program  refinements  are  suggested 
whenever  possible  to  aid  in  realizing  these  overall  objectives.  Finally, 
an  automatic  interrogation  "wrap-around"  program  is  described  for  uti¬ 
lizing  existing  missile  simulations  in  determining  the  required  incre¬ 
mental  equations  for  statistical  covariance. 


b.  Computer  Analysis  Objectives 

Basic  objectives  computer  analysis  tools  for  large- 
scale  systems  may  be  conveniently  5,rouped  under  the  interrelated  headings 
of  accuracy,  computational  speed,  computing  equipment  requirements,  and 
ease  of  implementation.  In  this  paragraph  the  statistical  covariance 
technique  will  be  examined  from  the  viewpoint  of  how  well  it  meets  these 
analysis  objectives. 

(1)  Accuracy .  The  primary  objective  of  any  computer 
analysis  technique  is  to  obtain  sufficient  accuracy  for  the  results  to 
be  trusted.  Without  some  degree  of  confidence  in  the  accuracy  of  a  given 
method,  one  cannot  proceed  with  assurance  to  the  more  complicated  design 
problem.  In  the  case  of  the  statistical  covariance  technique,  the  exact 
solution  is  obtained  for  linear  systems.  Furthermore,  a  very  accurate, 
though  not  exact,  solution  is  realized  for  mildly  nonlinear  systems  with 
low  noise  levels.  Only  for  highly  nonlinear  systems  and/or  high  noise 
levels  does  the  accuracy  become  unacceptable.  Even  in  those  cases,  the 
statistical  covariance  results  often  compare  favorably  with  the  results 
of  the  usual  25  to  50  Monte  Carlo  runs.  Moreover,  using  the  approach 
developed  in  Paragraph  3,  one  can  predict  to  some  extent  what  accuracy 
may  be  expected  by  statistical  covariance  on  a  given  problem.  Because 
the  statistical  covariance  technique  may  far  exceed  the  acceptable  level 
of  accuracy,  it  is  often  possible  to  operate  at  a  reduced  accuracy  to  aid 


in  achieving  one  or  more  of  the  other,  more  elusive,  interrelated  goals, 
such  as  computational  speed  or  core  requirements  of  the  digital  computer. 
For  example,  if  an  error  of  only  1  percent  is  being  achieved  by  sta¬ 
tistical  covariance,  it  might  be  advantageous  to  use  either  a  larger 
step  size  or  a  faster,  but  less  accurate,  integration  method  to  speed  up 
the  calculation  while  giving  an  error  of  5  percent. 

(2)  Computational  Speed.  The  second  objective  of  a 
computer  analysis  software  package  is  to  achieve  results  at  a  satis¬ 
factory  computational  speed.  As  mentioned  previously,  speed  is  closely 
related  to  accuracy  and  may  be  improved  if  less  accuracy  is  permitted. 

In  addition  to  increasing  speed  by  using  a  larger  step  size  or  a  faster 
integration  method,  e.g.,  Euler's  method,  there  are  certain  simplifying 
approximations  which  tend  to  speed  up  the  basic  statistical  covariance 
algorithm.  One  of  these  is  the  use  of  constant  coefficients  in  place 

of  slowly  time-varying  coefficients  in  the  incremental  equations.  These 
coefficients  vary  slowly  if  the  nominal  (noise-free)  trajectory  is  rela¬ 
tively  constant  over  a  time  span.  Moreover,  when  particular  coefficients, 
even  if  variable,  are  near  zero,  neglecting  them  entirely  can  often 
increase  computational  speed  significantly.  This  approximations  will  be 
utilized  in  the  automatic  sensitivity  program. 

(3)  Computing  Equipment  Requirements.  The  third  objec¬ 
tive  of  computer-oriented  tools  is  that  they  be  amenable  for  operating 
on  a  small  or  medium-sized  digital  computer,  i.e.,  the  limited  equipment 
aspect  of  all  engineering  problems  must  be  considered  in  introducing  a 
new  computer  analysis  tool.  One  consideration  is  word-length  of  the 
computer.  The  statistical  covariance  technique  can  be  easily  handled  on 
the  32-bit  Sigma  5.  Work-length  problems  would  arise  in  this  connection 
only  if  the  technique  were  being  used  on  much  smaller  bit  machines, 
e.g.,  minicomputers.  A  second,  much  more  important  consideration  is  the 
amount  of  core  required  for  processing  data  in  the  statistical  covariance 
technique.  A  sequential  algorithm  is  described  in  Paragraph  5  to  aid 

in  reducing  excessive  core  requirements.  A  rapid  access  device  (RAD) 
may  be  used  as  an  extended  core  for  overlaying  problems.  An  RAD  makes 
available  1.5  million  words  of  core  in  addition  to  the  32  K  words  of  the 
main  memory.  It  should  be  noted  that  using  the  RAD  decreases  the  com¬ 
putational  speed  of  a  given  algorithm.  In  particular,  there  is  a  delay 
caused  by  an  average  latency  time  of  17  msec  and  a  data  transfer  from 
the  RAD  to  the  main  memory  (50  msec  required  to  transfer  1500  words). 
Therefore,  recognizing  that  only  limited  equipment  is  available  is 
Important  in  the  design  of  computer-oriented  algorithms.  The  statisti¬ 
cal  covariance  technique  can  be  tailored  for  operating  on  limited  com¬ 
puting  equipment. 

(4)  Ease  of  Implementation.  A  fourth  objective  of  com¬ 
puter  analysis  techniques  is  that  the  algorithms  be  easily  applied  from 
the  user's  viewpoint  to  a  given  problem.  Special  programing  or  the 
need  for  precalculated  data  detracts  from  any  computer  method.  For  the 
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statistical  covariance  technique,  the  structural  form  of  the  coefficients 
for  the  incremental  (or  sensitivity)  equations  may  be  determined  in 
advance  by  taking  partial  derivatives  as  needed.  An  alternate  approach 
based  on  automatic  machine  computation  will  follow. 


c.  An  Automatic  Sensitivity  Computer  Package 

The  three  basic  parts  of  the  statistical  covariance  package 
are  the  determination  of  the  incremental  equation  coefficients,  the  inte¬ 
gration  of  the  covariance  equations,  and  the  integration  of  the  nominal 
system  equations.  The  purpose  of  this  section  is  to  describe  how  the 
incremental  equation  coefficients  may  be  automatically  computed  from  an 
existing  simulation  program  of  the  noise-free  system. 

Recall  that  for  the  nonlinear  system  given  by  Equation  (2),  i.e., 

x  =  f(x,  w,  t) 

the  incremental  equation  coefficients  are  given  by  Equation  (23),  i.e.. 


A(t) 
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w 


(t) 
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X  =  £ N  (t> 
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The  partial  derivatives  in  Equation  (23)  may  be  computed  numerically  as 

f(i„  +  y  c)  ■  f(v  y  c) 

A(t)  =  - 

f 

B(t)  =  - 

Aw 


Ax 


(V  +  '  f(SN,  y  t) 
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where  t^  represents  the  mean  value  of  the  Input  disturbance  vector  w, 

which  previously  has  been  assumed  to  be  zero.  The  notations  Ax  and  Aw 
in  Equation  (39)  represent  small  perturbations  about  the  nominal  tra¬ 
jectory  Xjj(t).  Because  P(t)  is  being  computed  within  the  same  time 

advance  loop.  A*  and  Aw  may  be  selected  as  one-tenth  of  the  standard 
deviation  of  x(t)  and  w(t),  respectively.  Whenever  any  main  diagonal 

element  of  P(t)  is  near  zero,  a  lower  limit  on  Ax,  such  as  10"^,  should 
be  utilized. 


This  numerical  evaluation  of  partial  derivatives  may  also  be  applied 
for  estimating  the  accuracy  expected  in  the  analysis  of  nonlinear  systems. 
Equations  may  be  formed  similar  to  Equation  (39)  by  using  AX  and  Aw  as 
one-sigma  values,  rather  than  one-tenth  of  a  sigma  in  each  case,  i.e.. 


— +  V  V  *)  ‘  1(~N>  Y  fc)  "  A ax  +  Il(*N’  V  fc) 
\  +  V  *)  '  Y  *)  "  B(t)  +  —2 N ’  V  fc) 


(40) 


The  vectors  g^  and  g^  in  Equation  (40)  represent  the  higher  order  terms 

in  the  Taylor  series  expansion  of  f  about  (xN>  ,  which  in  this  case 

has  only  been  approximated  numerically.  A  vector  k  may  be  computed 
whose  components  represent  the  percent  of  the  incremental  equations 

caused  by  the  higher  order  terms,  where 


— li^— N’  V  C)  +  -2i(~N’  r*w*  C) 

ki  A(t)  ax  +  B(t )  owj 


100  percent 


These  components  of  k  enable  the  results  of  Paragraph  3  to  be  used  to 
obtain  an  estimate  of  the  accuracy  to  be  expected  from  the  statistical 
covariance  technique. 

Although  these  numerical  procedures  for  determining  A(t)  and  B(t) 
are  yet  to  be  verified  in  actual  nonlinear  system  simulations,  they 
represent  a  realistic  approach  toward  achieving  a  self-contained  computer 
software  package.  The  nominal  system  equations  of  existing  simulations 
can  be  placed  inside  the  complete  package  which  will  interrogate  the 
system  equations  as  required  with  small  perturbations  to  form  A(t)  and 
B(t) .  It  will  not  even  be  necessary  to  know  what  the  system  equations 
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are;  only  their  number  and  the  number  and  location  of  noise  inputs  are 
required.  Such  a  comprehensive  computer  software  package  for  statisti¬ 
cal  covariance  would  be  extremely  useful  for  a  wide  variety  of  missile 
systems  analysis  problems. 

d.  Summary 

Program  refinements  for  statistical  covariance  have  been 
described  from  the  viewpoint  of  their  impact  on  basic  computer  analysis 
objectives:  accuracy,  computational  speed,  equipment  required,  and  ease 
of  implementation.  Numerical  procedures  have  been  developed  for  deter¬ 
mining  the  coefficients  of  the  incremental  equations.  Furthermore,  an 
estimate  of  the  accuracy  of  the  statistical  covariance  technique  may  be 
obtained.  Finally,  the  philosophy  of  a  self-contained  computer  software 
package  has  been  advanced  to  virtually  eliminate  the  majority  of  imple¬ 
mentation  problems. 


5.  Sequential  Calculations  of  the  Statistical  Covariance  Equations 
a.  Introduction 

In  large-scale  air  defense  missile  systems  it  is  con¬ 
venient  to  be  able  to  perform  computations  sequentially  to  avoid  exces¬ 
sive  digital  computer  storage  requirements.  Hesarovic  [35,  36]  proposed 
the  concept  of  multilevel  systems,  whereby  complex  systems  are  partitioned 
into  simpler  subsystems  for  analysis  and  design  purposes.  Mesarovic 
used  a  hierarchy  of  system  models  to  form  a  stratified  description  of 
complex  systems.  Lefkowitz  [37]  described  how  the  multilevel  hierarchy 
approach  has  been  used  to  solve  particular  industrial  problems.  Noton 
[38]  has  applied  multilevel  systems  theory  to  derive  a  coordination  algo¬ 
rithm  for  a  number  of  subsystem  Kalman  estimators. 

A  sequential  algorithm,  related  to  the  multilevel  system  concept, 
is  developed  in  this  paragraph  for  applying  the  statistical  covariance 
technique  to  large-scale  missile  systems.  As  shown  in  Figure  23,  the 
overall  system  is  segmented  into  several  subsystems  interconnected  by 

feed-forward  and  feedback  paths.  It  is  believed  that  large-scale  systems  j 

can  be  handled  more  efficiently  by  the  sequential  algorithm  than  by  the 

basic  covariance  program.  j 


b.  The  Sequential  Approach 

Consider  again  the  basic  statistical  covariance  equation 
developed  in  Equation  (17)  of  Paragraph  3  and  repeated  here  for  convenience. 

P  =  AP  +  pTAT  +  BQBT 
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Figure  23.  The  Basic  Large-Scale  System  Showing  Individual 
Subsystem  Connections 


Equation  (17)  describes  the  evolution  of  the  state  covariance  matrix 
P(t)  as  time  increases  from  tQ.  Let  the  large-scale  system  be  segmented 

into  several  subsystems  as  shown  in  Figure  23.  In  the  subsystem  context, 
the  matrices  A,  B,  P,  and  Q  may  be  partitioned  as 
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following  two  assumptions  are  made  at  this 

point  : 

1) 

The  matrix  P  is  symmetric. 

i .e . ,  P  = 

T 

P  and 

in  particular, 

P.  .  = 
ij 

T 

P.  . 

Ji 

2)  Disturbances  are  uncorrelated  with  each  other  and  each  enters 
only  the  single  designated  subsystems  as  shown  in  Figure  23.  This  means 
that  B  and  Q  are  diagonal. 

Therefore,  Equation  (17)  may  be  expressed  as 


P. .  =  A. . P. ,  +  A. „P„ .  +  A.-P-.  +  .  .  .  +  A.  P  . 
ij  ll  lj  i2  2j  i3  3j  lN  Nj 


+  AT  +  P^  AT  + 

+  Pu  AjX  +  P2.  Aj2  +  P3.  Aj3  + 


T  T 

•  •  +  Pm-  A.m 
Ni  jN 


+  B..  Q, .  B,  .  5.  . 
li  ii  ii  ij 


(41) 


where  is  zero  if  i  *  j  and  unity  if  i  =  j.  Equation  (41)  may  be 
conveniently  grouped  as 
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i-1 


Aik  Pkj  +  Aii  Pij 


k=l 


N 

+  I  Atk  pkl 

k=i+l 


j-1  N 

Z  ki  Ajk  +  *ji  Ajj  +  Z,  ‘ki  "jk  ‘  “ii  ^ii  “ii  uij 
k=l  k=j+l 


^  T  T  T 

2,  +  Q-  B*<  ** 


for  i  =  1  .  .  .  ,  N  and  j  =  1  .  .  .  ,  N 


(42) 


T  T 

The  matrices  and  P  in  the  second  line  of  Equation  (42)  may  be 
replaced  by  P^  and  P^  ,  respectively,  because  P  is  symmetric.  The 
second  and  fifth  terms  in  Equation  (42)  are  the  only  ones  involving  P 


ij 


Moreover,  the  first  and  fourth  terms  have  as  coefficient  matrices  entries 
from  the  lower  left  of  the  main  diagonal  of  the  system  matrix  A.  These 
terms  represent  feed-forward  paths.  Elements  from  the  upper  right  of 
the  main  diagonal  of  A  appear  in  the  third  and  sixth  terms  in  Equation 
(42)  and  represent  feedback  paths.  The  seventh,  or  last,  term  represents 
input  noise  data.  Rewriting  Equation  (42)  and  summarizing  these  obser¬ 
vations,  yields 


Self- Interact ion 


Feed -Forward  Paths 


N 

I 

[_k=i+l 


Aik  Pkj  + 


N 

I 

k=j+l 


T 

P  A 
ik  Ajk 


+  Bii  Qii  Bii  6ij 


(43) 


Feedback  Paths 


Noisy  Inputs 


What  is  desired  is  to  apply  Equation  (43)  sequentially  to  determine  P 

for  all  i  and  all  j.  It  is  particularly  important  only  to  know  P  ,  but 

it  may  be  shown  that  calculation  for  all  i  and  j  is  necessary  to  com¬ 
pletely  determine  P^. 
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For  example,  consider  the  problem  of  determining  P^  for  the  case 
where  no  feedback  terms  are  present.  Using  Equation  (43)  yields 

*44  =  A44  P44  +  P44  A44  +  A43  P34  +  P34  A43  +  B44  Q44  ®44 


P34  =  A33  P34  +  P34  A44  +  A32  P24  +  P33  A43 


P24  A22  P24  +  P24  A44  +  A21  P14  +  P23  A43 


P14  “  A11  P14  +  P14  A44  +  P13  A43 


(44) 


In  interpreting  Equation  (44),  note  that  the  four  matrix  equations  must 
be  applied  sequentially  in  reverse  order,  beginning  with  the  last.  The 
matrices  P^>  P23>  and  P33  are  known  calculations  for  the  previous 

subsystem,  i.e.,  No.  3.  It  is  assumed  that  subsystem  No.  4  has  inputs 
only  from  subsystem  No.  3  and  from  external  sources. 

A  flow  chart  describing  the  sequential  method  for  statistical 
covariance  calculations  is  provided  in  Figure  24.  Additional  investi¬ 
gation  is  required  to  develop  a  workable  computer  software  package  for 
use  in  large-scale  missile  system  analysis. 


c .  Summary 

This  paragraph  has  described  a  sequential  algorithm  as 
a  means  of  handling  large-scale  analysis  problems  more  efficiently  than 
the  basic  statistical  covariance  program.  The  new  algorithm  categorizes 
terms  on  a  subsystem  basis  as  self- interacting,  feed- forward  paths, 
feedback  paths,  and  noisy  inputs.  Further  work  is  needed  to  develop 
the  sequential  algorithm  into  a  useable  computer  software  package. 

6.  Conclusions  and  Recommendations 
a.  Conclusions 

The  primary  conclusion  of  this  report  is  that  the  statis¬ 
tical  covariance  technique  yields  significantly  faster,  more  accurate, 
results  than  the  traditional  Monte  Carlo  approach  for  a  wide  range  of 
missile  systems  analysis  error  propagation  problems.  The  statistical 
covariance  algorithm,  which  gives  exact  results  for  linear,  time-varying 


systems,  may  be  applied  to  nonlinear  systems  by  using  linearized  incre¬ 
mental  equations  about  the  noise-free  solution.  It  was  stated  in  Para¬ 
graph  3  that  the  accuracy  of  the  statistical  covariance  technique  for 
nonlinear  systems  depends  entirely  upon  the  relative  accuracy  of  this 
linearizing  approximation.  A  parameter  (a)  of  a  second-order  nonlinear 
example  was  varied  to  demonstrate  that,  as  a  system  becomes  highly  non¬ 
linear,  a  larger  error  is  obtained.  Moreover,  as  the  input  noise  level 
( Q )  was  increased  to  very  high  values,  the  error  in  the  statistical 


covariance  results  was  increased.  It  may  be  concluded  that,  while  there 
admittedly  are  extreme  conditions  for  which  the  statistical  covariance 
algorithm  yields  unacceptable  errors,  the  technique  works  satisfactorily 
for  the  conditions  occurring  in  most  missile  system  analysis  problems. 


A  special  interrogation  program  has  been  suggested  for  automatically 
computing  the  incremental  equation  coefficients  by  using  existing  missile 
system  simulations.  This  "wrap-around"  program  interrogates  an  existing 
simulation  program  to  determine  locally  linear  sensitivity  coefficients- 
Therefore,  it  is  unnecessary  to  reprogram  a  complex  missile  system  simu¬ 
lation  to  apply  the  statistical  covariance  technique.  Furthermore,  the 
automatic  sensitivity  program  is  capable  of  computing  information 
required  for  estimating  the  accuracy  of  the  statistical  covariance  results 
and  the  number  of  Monte  Carlo  runs  needed  for  comparable  accuracy. 


Procedures  for  transforming  the  statistical  covariance  results  to 
determine  miss  distances  were  also  described  in  Paragraph  3.  A  short 
integration  program  is  required  for  handling  the  results  of  the  trans¬ 
formation  for  a  complete  missile  systems  analysis  package. 


A  sequential  algorithm  was  developed  in  Paragraph  o  to  aid  in 
reducing  the  amount  of  required  computer  storage.  The  new  algorithm 
treats  the  overall  problem  on  a  subsystem  basis  and  performs  statistical 
covariance  calculations  serially  from  subsystem  to  subsystem.  An  addi¬ 
tional  approximation  on  sampling  in  feedback  loops  must  be  made  to  permit 
the  sequential  ordering  of  these  calculations. 


It  has  been  demonstrated  that  the  Monte  Carlo  approach  often  yields 
poor  results  in  error  propagation  studies.  In  Paragraph  2  effort  was 
devoted  to  obtain  a  correct  implementation  of  the  Monte  Carlo  method. 

In  particular,  several  pseudo-random  number  generators  were  tested  on  a 
second-order  linear  system.  It  was  noted  that,  even  for  the  best  of 
these  generators,  1000  Monte  Carlo  runs  only  yielded  an  average  error  of 
1.5  percent  for  the  output  variance.  Furthermore,  for  the  same  linear 
system  errors  of  30,  15,  and  10  percent  were  obtained  for  25,  50,  and 
100  runs,  respectively.  Similar  results  were  obtained  for  a  second- 
order  nonlinear  system.  The  conclusion  is  that  a  larger  number  of  Monte 
Carlo  runs  than  usually  performed  are  required  for  acceptable  accuracy. 
Rather  than  expending  excessive  computer  time  for  large  numbers  of  Monte 
Carlo  runs,  a  more  satisfactory  approach  is  to  utilize  the  faster  statis¬ 
tical  covariance  technique  in  a  totally  preprogrammed  computer  software 
package . 
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It  has  been  shown  that  the  statistical  covariance  technique  is 
applicable  to  a  wide  range  of  missile  system  analysis  problems.  In 
Paragraph  4, the  technique  was  evaluated  according  to  four  general  objec¬ 
tives  of  all  computer  analysis  tools:  accuracy,  computational  speed, 
equipment  requirements,  and  ease  of  implementation.  Extreme  accuracy 
was  demonstrated  in  Paragraphs  2  and  3  on  linear  and  nonlinear  sys¬ 
tems.  The  accuracy  prediction  scheme  may  become  important  for  highly 
nonlinear  systems  with  high  noise  levels.  The  computational  speed  was 
compared  favorably  with  Monte  Carlo  in  Paragraph  4.  The  sequential 
algorithm  of  Paragraph  5  was  developed  to  help  in  solving  the  computer 
storage  problem,  and  the  use  of  RAD  further  aids  in  making  equip¬ 
ment  requirements  less  critical.  Finally,  the  implementation  problem 
is  greatly  reduced  by  using  the  automatic  sensitivity  program  for  coeffi¬ 
cient  calculations  as  an  integral  part  of  the  computer  software  package. 


b.  Recommendations  for  Further  Work 

The  first  recommendation  is  that  a  total  computer  software 
package  for  statistical  covariance  be  developed  especially  tailored  for 
missile  systems  analysis  problems.  In  addition  to  the  basic  algorithm, 
this  package  should  include  the  automatic  sensitivity  program  for  cal¬ 
culating  incremental  equation  coefficients,  an  accuracy  prediction  pro¬ 
gram  associated  with  the  automatic  sensitivity  program,  and  a  miss  dis¬ 
tance  transformation  program.  A  comprehensive  self-contained  package 
containing  these  programs  would  require  only  minimal  information  from 
the  user  while  providing  the  flexibility  of  performing  difficult  missile 
system  analysis  production  runs. 

The  second  recommendation  is  that  the  sequential  version  of  the 
statistical  covariance  technique  be  programned  and  tested  as  a  means  of 
reducing  computer  storage  requirements.  Only  when  it  has  been  shown 
that  a  significant  savings  in  storage  is  realized  by  the  sequential 
algorithm  should  it  be  incorporated  into  the  total  computer  software 
package . 

Finally,  the  third  recommendation  is  that  improvements  be  made  in 
Monte  Carlo  simulations.  The  major  deficiency  at  present  is  the  need 
for  a  better  random  number  sequence  to  serve  as  system  inputs.  In 
addition,  more  work  on  the  discretization  problem  of  continuous  random 
signals  could  be  useful. 

The  development  of  these  improvements  in  statistical  techniques 
would  greatly  enhance  missile  system  analysis  capabilities.  Better 
analysis  techniques,  as  those  described  herein,  must  be  further  developed 
and  tested  in  practical  systems  operations. 
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Appendix  A.  COMPUTER  PROGRAMS  FOR  PSEUDO-RANDOM 
NUMBER  GENERATION 


This  appendix  provides  additional  computer  programs  and  output  data 
on  the  generation  and  testing  of  pseudo-random  numbers  from  multiplica¬ 
tive  formulas.  Figure  A-l  gives  the  output  of  the  computer  program  of  the 
20 

Brown,  (M  ■  2  )  generator  shown  in  Figure  1.  Corresponding  results 

31 

for  the  Brown  (M  =  2  )  generator  are  included  in  Figures  A-2  and  A-3. 

Figure  A-4  shows  the  computer  program  used  to  obtain  information  for  the 
20 

M  «  2  curves  in  Figures  2  and  3.  The  program  in  Figure  A-5  was  used 
for  the  Monte  Carlo  results  plotted  in  Figures  5  and  6.  These  detailed 
computer  programs  have  been  included  in  this  Appendix  for  use  in 
the  exact  reproduction  of  the  two  Brown  multiplicative  generators 
described  in  Paragraph  2  and  their  use  in  Monte  Carlo  simulations. 
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Figure  A-l.  The  First  30  Numbers  from  the  Brown  (M  =  2  )  Generator 

in  Figure  1 
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1  C  OPERATION  OF  A  MULTIPLICATIVE  DIGITAL  GENERATOR 

2  C 

3  C 

4  XMEAN=0 . 

5  SIG=1 - 

6  IX=31571 

7  DUM=0 . 1 

8  DO  99  K=l,30 

9  C 


10 

C 

MULTIPLICATIVE  GENERATOR 

11 

C 

WITH  M=2  TO  THE  3 1ST  POWER 

12 

C 

13 

IY=19971*IX 

14 

IX=IY 

15 

IF  (IX)  5,5,6 

16 

5  IY=IY+2 14 7483647+1 

17 

6  U=IY*0.4656612873E-9 

18 

C 

19 

c 

TRANSFORMATION  TO  NORMAL 

20 

c 

21 

Z=SQRT  (-2 . 0*AL0G (DUM) )*SIG 

22 

XNORM=Z  *COS (6 . 28318*U)+XMEAN 

23 

c 

24 

DUM=U 

25 

PRINT  77,K,IX,U,XNORM 

26 

77  FORMAT (6X, 14, 8X , 112 , 2F15 .6) 

27 

99  CONTINUE 

28 

STOP 

29 

END 

Figure  A-2.  The  Brown  Multiplicative  Pseudo-Random  Number  Generator 
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Figure  A -4.  A  Digital  Computer  Program  for  Testing  the  Normally 
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Distributed  Brown  (M  =  2  )  Generator 
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42  CONTINUE 

32  CONTINUE 
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00  62  NA»1.MTPT 
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Figure  A-5.  A  Monte  Carlo  Simulation  Program  for  the  System 


in  Equation 


(7)  Using  the  Brown  (M  =  2 


Generator 


58 


5? 

53 

54 

55 

56 

57 

58 

59 

50 

51 
5? 
63 
5* 

55 

56 

57 

58 

59 


1 

p 

3 

4 
K 

5 


T«H#X‘!A*1C» 

SOL (N4/KA ) *.08333333333  -0 • 5 *FXP ( «2« 0*T ) 

1  ♦0.5A6*6*667*EXP< - T.O#T)  - " . 25*EXP ( -4 • 0*T ) 

S(NA,k'A)*3(NA,<A}/XNIIJF 

PIF  (* A»KA  )«l''f'.«(S(NA.KA  )-35Lr'A/K.A  )  ) /SGL  ( *-)A  ,  K  A  ) 
PRINT  /<T.S(KIA#<A  )  t  S^l  f  NA.KA  )  ,TIF  (NA.XA  ) 

7  FORMAT ( 10X.F5, 2. 3F15»6> 
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STOP 
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Figure  A-5.  Concluded 
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Appendix  B.  THE  EXACT  CONVOLUTION  SOLUTION  FOR  THE  OUTPUT 
VARIANCE  IN  LINEAR  SYSTEMS 


Consider  a  linear  system  with  zero  input  for  t  <  0  and  with  a  sta¬ 
tionary  random  process  w(t)  as  input  for  t  5  0.  Let  w(t)  be  a  zero- 
mean  white  noise  process  with  variance  Q^.  Because  the  input  statistics 

are  different  for  the  two  ranges  of  t,  the  effective  input  must  be  con¬ 
sidered  to  be  nonstationary.  Papoulis  [39]  has  shown  that  the  expression 
for  the  output  autocorrelation  is  given  by 

*«(*!•  h)  ■  Ew(V  e2)  *  h(K)  *  h  ('*)  '  <B-1) 

where  *  represents  the  convolution  operation  and  h(t)  is  the  impulse 
response  of  the  linear  system,  i.e.,  the  inverse  Laplace  transform  of 
the  transfer  function  H(s). 


As  an  example,  consider  the  second-order  system  given  in  Equation  (7) 
of  Paragraph  1.  The  impulse  response  h(t)  is  obtained  by  using  Equation 
(8),  i.e., 


h(t)  =  L-1  H(s)  =  L"1 


(s  +  1)  (s  +  2) 


=  g-t  .  e'2t  t  >  o 


(B-2) 


Performing  the  first  convolution  in  Equation  (B-l)  yields 


’»('!'  C2)  *  h(tl)  ‘  V(‘l  -  C2)  *  (*'“  -  e'2t0 

for  tj^  ^  0,  t2  ^  0 

- /  v(‘i  -  *2  -  ’)  KT  -  •■*]  * 

o 

(  ^[e’  (tl  ’  t2)  -  e-2^1  ’  t2)]  i  h  -  fc2  2  0 

J  0  Otherwise 

Pteceiinj  page  blank  '  < 
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The  second  convolution  yields 


***(*1’  C2)  =  h)  *  h(tl)]  *  h(C2) 

=  f  Q^e'fc  “  t2  +  T)  -e'2^1  "  t2  +  T)]  •  [e-*  -  e’2*]  dT 
J0 

•■(‘i-  '2)  - 1.  °<C1  ‘ t2)  [1  (l  -  «'2t2)  -  j  (l  -  *‘3t2)] 

-  »„e-2(n  -  '2)  [I  (,  .  .-**)  -  i  (l  -  .-^2)] 

for  tj  i  t,  s  0  .  (B-4) 

Therefore,  the  output  variance  Q^(t)  is  given  by 

Qx(t)  =  Rxx(Cl’  *"2) 

.  n  T 1  1  -2t  ,  2  -3t  1  -4t 

Qx(t)  =  %\U  '  2  e  +  3  G  *  4  e  J 

for  t  -  0  .  (B-5) 

This  exact  solution  is  used  for  comparison  with  Monte  Carlo  simulations 
in  Figure  5 . 
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Appendix  C.  COMBINED  STATISTICAL  COVARIANCE  AND  MONTE  CARLO  I 

COMPUTER  PROGRAMS 


The  combined  statistical  covariance  and  Monte  Carlo  simulation  pro¬ 
gram  for  a  second-order  nonlinear  system  (Figure  C-l)  is  included  in  this 
appendix.  This  basic  program  was  modified  for  different  input  data  to 
produce  the  numerical  results  given  in  Figures  11  through  20. 
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C - X  "COMBINED  PHOORAH  FSB  STAT  IgTlCAITCOVARI  ANCE - - 
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Figure  C-l.  A  Combined  Program  for  Statistical  Covariance 
and  Monte  Carlo  for  a  Nonlinear  System 
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Figure  C-l.  Continued 
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